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Abstract. We propose a mixed finite element method for the motion of a 
strongly viscous, ideal, and isentropic gas. At the boundary we impose a 
Navier-slip condition such that the velocity equation can be posed in mixed 
form with the vorticity as an auxiliary variable. In this formulation we design 
a finite element method, where the velocity and vorticity is approximated with 
the div- and curl- conforming Ncdclec elements, respectively, of the first order 
and first kind. The mixed scheme is coupled to a standard piecewise constant 
upwind discontinuous Galcrkin discretization of the continuity equation. For 
the time discretization, implicit Euler time stepping is used. Our main result is 
that the numerical solution converges to a weak solution as the discretization 
parameters go to zero. The convergence analysis is inspired by the continuous 
analysis of Feireisl and Lions for the compressible Navier— Stokes equations. 
Tools used in the analysis include an equation for the effective viscous flux 
and various renormalizations of the density scheme. 
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1. Introduction 



Let ft C R , N = 2, 3, be an open, convex, polygonal domain with Lipschitz 
boundary dfl and let T > be a final time. We consider the flow of an ideal 
iscntropic viscous gas governed by the Stokes approximation equations 



Here, the unknowns are the density g = g(t, x) > and velocity u = u(t, x) £ R N . 
The operators V x and div^ are respectively the spatial gradient and divergence 
operators, and A = div x V x is the Laplace operator. The viscosity coefficients /z, 
A are assumed to be constant and to satisfy fj, > 0, NX + 2/i > 0. 

The pressure is given by Boyle's law which in the iscntropic regime takes the 
formp(g) = ag 1 , where a > is constant. In real applications the value of 7 ranges 
from a maximum of | for monoatomic gases, to values close to one for polyatomic 
gases at high temperatures. In this paper, we will for purely technical reasons be 
forced to require 7 > ^. 

From the point of view of applications, the model (f.l)-(f.2) can be justified 
for flows at very low Reynolds numbers so that the effects of convection may be 
neglected. It is also on the same form as various shallow water models [11]. From 
a mathematical perspective, the system (1.1)-(1.2) is a model problem containing 
some, but not all, of the difficulties associated with compressible fluid dynamics. 

Mathematical analysis concerning the well-poscdness of the system (1.1)— (1.2) 
seems to originate with the papers [12, 14] by Kazhikov and collaborators. Several 
other contributions on the existence and long term stability exist, also in the context 
of similar shallow water models. However, for our purpose here, the most relevant 
study is that of Lions [11] in which the global existence of (weak) solutions and 
some higher regularity results are established. 

In this paper we impose the following boundary conditions: 



where v is the unit outward normal on dfl and curLj is the curl operator. Here, 
in 2D, curia; denotes the rotation operator taking vectors into scalars. The first 
condition is a natural condition of impermeability type on the normal velocity. 
The second condition is in the literature commonly referred to as the Navier-slip 
condition. While these boundary conditions are not motivated by physics, they 
are widely used in numerical methods. In particular, in the context of geophysical 
flows they are often preferred over classical Dirichlet conditions since the latter 
necessitates expensive calculations of boundary layers. Of more importance to this 
paper, the boundary conditions (1.3)- (1.4) allow us to pose the system (1.1)— (1.2) 
in mixed form with curL u as an auxiliary variable. This fact will play a crucial 
role in the upcoming analysis. 

While many numerical methods appropriate for the Stokes approximation and 
Navier-Stokes equations have been proposed, the convergence properties of these 
methods arc mostly unsettled. In fact, it is not clear whether or not any of these 
methods, in more than one dimension, converge to a (weak) solution as discretiza- 
tion parameters tend to zero. In one dimension, there are some available results 
due to D. Hoff and his collaborators. However, these results apply to an ideal gas 
in Lagrangian coordinates and with initial data of bounded variation. In more 
than one dimension, there are some recent results for simplified models. In the pa- 
pers [6,7], a convergent finite element method for a stationary compressible Stokes 



d t g + div x (gu) = 0, in (0, T) x il, 
d t u — fxAu — AVz divz u + V x p{g) = 0, in (0, T) x fl. 



(1.1) 
(1.2) 



u-u = 0, on(0,T)x<9Q, 
curlz u x v = 0, on (0, T) x dfl, 



(1.3) 
(1.4) 
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system is proposed and analyzed. The system considered there are similar to (1.1)- 
(1.2) but without temporal dependence. In [9], we established convergence of a 
finite element method for a semi-stationary version of (1.1)— (1.2) (dtu = 0) and 
homogenous Dirichlet boundary conditions. This paper can be seen as a continua- 
tion of the recent study [8] in which a convergent numerical method for the same 
semi-stationary system ((1.1)— (1.2) with dtu = 0) with boundary conditions (1.3)- 
(1.4) was established. The main novelty of this paper is consequently the addition 
of the time derivative term dtu in the velocity equation (1.2). 

Let us now discuss our choice of numerical method for the Stokes approximation 
equations. For the time discretization, we will use implicit time stepping in both 
equations. To approximate the continuity equation (1.1) we will use a standard 
piecewise constant upwind discontinuous Galerkin method. To approximate the 
velocity, we will use a mixed finite element method with the Nedelec's spaces of the 
first order and first kind. The mixed formulation is motivated by introducing the 
vorticity w — curb; u as an auxiliary unknown and recasting the velocity equation 
(1.2) in the form: 

diu + ficm\ x w — (A + p)V x divz u + V x p(g) = 0, (1-5) 

where the identity — A = curls curl x — V x div^ is used. This leads to a natural mixed 
formulation in which the requirement w = curl^ u plays the role of a lagrangian 
multiplier. 

Denote by Wq 1V ' 2 (£Y) the vector fields it on O for which div x u G L 2 and u ■ 
v\dn = 0, and by W^ url ' 2 (f2) the vector fields w on £1 for which curl x w £ L 2 and 
w x v\gn = 0. We choose corresponding finite element spaces Vh C Wq lv,2 (f2) and 
Wh C W =url,2 (ri) based on Nedelec's elements of the first order and first kind [13]. 
The mixed finite element methods seeks, for each time step k — 1, . . . , M , functions 
(«>*,«*) G Wh x Vh such that 



(1.6) 



&t (u k h ) v h + /x curl, w h h v h + [(A + /a) div x u\ - p{g k h ) div x v h ] dx = 0, 

/ w*rjh - u\ curia, t] h dx = 0, 
Jq 

for all (v h ,r) h ) G W h x V h , where q\ is given and d£ (u%) = (Ai) -1 [w£ - u£ -1 ] 
denotes implicit time stepping. Note that the boundary conditions (1.3)-(1.4) are 
mandatory to obtain this formulation. 

Our main result is that {(wh, Uh, Qh)}h>a converges to a weak solution of the 
Stokes approximation equations, at least along a subsequence. The major difficulty 
is to obtain strong compactness of the density approximation {gh}h>o which is 
needed in order to pass to the limit in the nonlinear pressure function. Since the 
density approximations arc only bounded in L°°(0, T; L 7 (f2)) this is intricate. At 
the heart of the convergence analysis lies the effective viscous flux P e s(Qh-,Uh) = 
p{Qh) — (A + fi) div x Uh- In particular, strong convergence of the density approxi- 
mation follows from the property: 

lim J J ipP cS (Qh,Uh)gh dxdt = J J \pP cS (g,u)g dxdt, (1.7) 

for all ip G C£°(0,T). It is in the process of obtaining (1.7) that the carefully 
selected finite element spaces and mixed form prove useful. Specifically, we obtain 
(1.7) by setting Vh = II^VajA £>h in (1.6), where 11^ is the canonical interpolation 
operator into Vh ■ This test function satisfies div x Vh = Qh and is almost orthogonal 
to curls. The main difficult in obtaining (1.7) is to treat the time derivative term, 
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which, with Vh as described above, is of the form 

J /^(w^n^VxA- 1 ^] dxdt = J | af wVxA- 1 ^] dxdt + o(h) 

= J J A-Mdiv^uhl^tefc) dxdt + 0{h). 

Using the continuity scheme, the last term the last term can be shown to converge. 
The property (1.7) then follows. Our analysis resembles that of Lions and Feireisl 
for the compressible Navicr-Stokes equations. 

As part of the analysis, we will need that the discrete velocity Uh converges 
strongly to a function u. This is not immediate since the approximation space 
Vh is only div^ conforming. To obtain strong convergence, we utilize the discrete 
Hodge decomposition Vh = curl x Wh + V h °' ± satisfied by the chosen Nedelec spaces. 
When writing Uh = curLj Ch + Zh, it can be seen that Ch does not depend on the 
density Qh and as a consequence converges strongly. The remaining term Zh is 
then weakly discrete curl free with bounded divergence and an estimate from the 
previous paper [8] yields 

\\z h (t,x)-Zh(t,x-Z)\\ L 2( 0tT . L 2 {Q - )) ^>0, as |£| ->0, (1.8) 

uniformly in h. From the velocity scheme, we deduce a weak time-continuity esti- 
mate of the form 

t h (z h )GL 1 (O,T;W- 1 ' 1 (n)) J (1.9) 

independently of h. The two estimates (1.8) and (1.9) tells us that Zh satisfies 
the hypotheses of an Aubin Lions type lemma (see Lemma 2.3 below for details). 
Strong convergence of Zh follows from this lemma. 

The paper is organized as follows: In Section 2, we introduce notation and list 
some basic results needed for the later analysis. Moreover, we recall the usual 
notion of weak solution and introduce a mixed weak formulation of the velocity 
equation. Finally, we introduce the finite element spaces and review some of their 
basic properties. In Section 3, we present the numerical method and state our main 
convergence result. Section 4 is devoted to deriving basic estimates. In Section 5, 
we establish higher intcgrability of the density. Finally, in Section 6, we prove the 
main convergence result stated in Section 3. The proof is divided into several steps 
(subsections), including convergence of the continuity scheme, weak continuity of 
the discrete viscous flux, strong convergence of the density approximations, and 
convergence of the velocity scheme. 



2. Preliminary material 

We will write W m ' p (Vl) for the Sobolev space of functions with derivatives of 
all orders up to m belonging to the space L p (Vl). To distinguish between scalar 
and vector functions, we will write vector functions with a bold face. Similarly, a 
functions space written in bold face denotes the vector analog of the corresponding 
scalar space. 

We make frequent use of the divergence and curl operators and denote these 
by div-E and curl x , respectively. In the 2D case, we will denote both the rotation 
operator taking scalars into vectors and the curl operator taking vectors into scalars 
by curia. 

We will make use of the spaces 

L 2 (n) = (</>€ L 2 (n) : J <f> dx = j , 
W div ' 2 (fl) = {ve L 2 (n) : &iv x v G L 2 (fl)} , 
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W cuxl ' 2 (Q) = {veL 2 (Q): curl x v G L 2 (Q)} , 

where v denotes the outward pointing unit normal vector on <9f2. If v € W dlv,2 (£l) 
satisfies v ■ v\qq_ = 0, we write v e W lv ' 2 (fl). Similarly, v G WQ Ur,2 (f2) means 
v G W curl,2 (0) and v x v\dn = 0. In two dimensions, to is a scalar function and 
the space Wg Ur ,2 {fl) is to be understood as Wq ,2 {Vl). To define weak solutions, we 
shall use the space 

W(fl) = {v G L 2 (n) : div x v G L 2 (Q),cut1 x v G L 2 (ft),u • i/| an = 0} , 

which coincides with W dIV ' 2 (ft) n W curM (O). The space W(fi) is equipped with 
the norm \\v\\ 2 w = \\v\\ 2 L2{n} + \\div x v\\ 2 L2{n) + ||curl x v\\ 2 L2(n) . It is known that ||-|| w 
is equivalent to the W 1 - 2 norm on the space {v G W 1,2 (f2) : v ■ v\en = 0}, see, e.g., 
[11]. 

For the convenience of the reader we list some basic functional analysis results 
to be utilized (often without mentioning) in the subsequent arguments (for proofs, 
see, e.g., [5]). Throughout the paper we use overbars to denote weak limits. 

Lemma 2.1. Let O be a bounded open subset ofM. M with M > 1. Suppose g: R — * 
(—00,00] is a lower semicontinuous convex function and {v n } n>1 is a sequence of 
functions on O for which v n — v in i 1 (0), g(v n ) € i 1 (0) for each n, g(v n ) 
g(v) in L l (0). Then g(v) < g(v) a.e. on O, g{v) G L l (0), and J Q g(v) dy < 
liminfn^oo J Q g(v n ) dy. If in addition, g is strictly convex on an open interval 
(a, b) C K and g(v) = g(v) a.e. on O, then, passing to a subsequence if necessary, 
v n (y) -> v(y) for a.e. y E {y G O \ v(y) G (a, b)}. 

Let X be a Banach space and denote by X* its dual. The space X* equipped 
with the weak-* topology is denoted by X^, oak , while X equipped with the weak 
topology is denoted by X wea k- By the Banach- Alaoglu theorem, a bounded ball in 
X* is (?{X* , X)-compact. If X separable, then the weak-* topology is metrizable 
on bounded sets in X*, and thus one can consider the metric space C ([0, T]',X^ cak ) 
of functions v : [0, T] — ► X* that are continuous with respect to the weak topology. 
We have v n -> v in C([0,T];X* eak ) if (v n (t),4>)x*,x -> (v{t),4>) x *,X uniformly 
with respect to t, for any cf> G X. The following lemma is a consequence of the 
Arzela-Ascoli theorem: 

Lemma 2.2. Let X be a separable Banach space, and suppose v n : [0, T] — > X* , 
n = 1, 2, . . . , is a sequence for which \\v n \\ Loo n xyx*) — ^> f or some constant C 
independent of n. Suppose the sequence [0,T] 3 t (v n (t),$)x*,x, n = 1,2,..., 
is equi- continuous for every $ that belongs to a dense subset of X. Then v n belongs 
to C ([0, T]; X£ eak ) for every n, and there exists a function v G C ([0, T}; A* cak ) 
such that along a subsequence as n — > 00 there holds v n — > v in C ([0, T]; A* eak ). 

In what follows, we will often obtain a priori estimates for a sequence {f n }„>i 
that we write as a v n €b A" for some functional space X. What this really means 
is that we have a bound on ||fn|lx that is independent of n. 

The following discrete version of a lemma due to Lions [11, Lemma 5.1] will 
prove useful in the convergence analysis. A proof of this lemma can be found in [9]. 

Lemma 2.3. Given T > and a small number h > 0, write (0,T] = U^L^tk-i, tk\ 
with tk = hk and Mh = T. Let {/h}^>o, {ffh}/^o ^ e ^ wo sequences such that: 

(1) the mappings t — > gh(t,x) and t — * fh(t,x) are constant on each interval 
(ifc-i.tfcj, * = !,..., AT. 



() 



KENNETH H. KARLSEN AND TRYGVE K. KARPER 



(2) {fh}h>o and {gh}h>o converge weakly to f and g in L Pl (0, T; L qi (fi)) ant 
L P2 (0, T; L q2 (Q)), respectively, where 1 <p\,q\ < oo and 



1 1 _ 1 1 

Pi P2 qi q_2 

(3) the discrete time derivative satisfies 

9h (t,x)- 9h (t-h,x) &bL x^ T . w -i,x m 
h 

(4) \\fh{t,x) - fh{t,x- 0\\lv2(p,T;L«2(Q)) -* as |£| -» 0, uniformly in h. 
Then gnfh ~^ gf * n the sense of distributions on (0, T) x Q. 

2.1. Weak and renormalized solutions. 

Definition 2.4 (Weak solutions). We say that a pair (g, u) of functions constitutes 
a weak solution of the Stokes approximation equations (1.1)— (1.2) with initial data 

( 5 ° )tl °)eF(n)xl 2 (fl), 7 >^, 

and Navier-slip type boundary conditions (1.3)-(1.4), provided the following con- 
ditions hold: 

(1) ( ft «) G L°°(0, T; LT(n)) x L 2 (0, T; W) n L°°(0, T; L 2 (ft)); _ 

(2) £> t + div^pw) = in the weak sense, i.e, V</> G C°°([0,T) x Q), 

f f g{cj)t+u\7 x (j)) dxdt+ f eV|t=o dx = 0; (2.1) 
Jo Jn Jn 

(3) « t - ij.Au - XW X div x u + W x p(g) = weakly, i.e, V0 G C°°([0, T) x IT) for 
which (f> ■ v = on (0, T) x 90, 

r-T 

—u<pt + /-t curl,,; it cur^ <A 

(2-2) 

+ [((j, + A) diva; u — p(g)] div x <p dxdt = / u°<p\ t= o dx. 

Jn 

For the convergence analysis we shall also need the DiPerna-Lions concept of 
renormalized solutions of the continuity equation. 

Definition 2.5 (Renormalized solutions). Given u G L 2 (Q, T; VV(fi)), we say that 
p G L°°(0, T; L 7 (f2)) is a renormalized solution of (1.1) provided 

B(g) t + div x (B(g)u) + b(g) div x u = in the weak sense on [0, T) x il, 

for any B G C[0, oo) n C^oo) with B(0) = and 6(g) := ^'(g) - 5(g). 

We shall need the following lemma. A proof can be found in [8]. 

Lemma 2.6. Suppose (g,u) is a weak solution according to Definition 2.4- If 
g G L 2 ((0,T) x fi)), then g is a renormalized solution according to Definition 2.5. 

2.2. A mixed formulation. In view of the Navier-slip boundary condition (1.4) 
the velocity equation (1.2) admits the following mixed weak formulation, which we 
will use to design a mixed finite clement method: Determine functions 



(w,u) G L 2 (Q,T; Wr l2 m x L 2 (0, T; W div ' 2 (S!)) 



such that 

i-T 



I I —uvt + /x curl x wv + [(/x + A) div^ u — p(g)] div x v dxdt = / u°v\ t=0 dx, 
Jo Jn Jo. 

I I wr) — curlj; rju dxdt = 0, 



(2.3) 
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for all {t],v) G L 2 (0,T;Wa' 2 (n)) x L 2 (0, T; W (O)) n W l2 (Q, T; L 2 (Q)). 

Note that if (w, u, g) is a triple satisfying the mixed formulation (2.3) then the 
pair (it, g) satisfies the weak formulation (2.4) [8]. 

2.3. Finite Element spaces and basic results. Throughout this paper, {Eh}h 
denotes a shape regular family of tetrahedral meshes of f2, where h is the maximal 
diameter. By shape regular we mean that there exists a constant k > such that 
every E G Eh contains a ball of radius \e > — , where He is the diameter of E. 
For each fixed h > 0, we let I\ denote the set of faces in Eh and Vh the set of 
edges. In two dimensions, Th is the set of edges and Vh the set of vertices. We will 
use Fj(E) to denote the space of vector polynomials on E with I components and 
maximal order k. 

To approximate the vorticity w, we will use the curl-conforming Nedelec space 
of the first order and kind [13]: 

W h ((l) = [w G W curl ' 2 (fi) : w\ E G W(£), V£ £ E h) 

/• 1 (2-4) 

y [to • t] e dS(x) = 0, Ve G V,, 



where f is the unit tangential along the edge e, [-] e is the jump over the edge e, and 
W{E) 



F\(E), N = 2, 

{io G P?(£) : V x it> + V a ii; T = 0} , AT = 3. 



In two dimensions, the continuity requirement J \w ■ i] dS(x) = in (2.4) is to 
be understood as continuity at vertices. For the velocity it, we will use the div- 
conforming Nedelec space of the first order and kind [13]: 

V h (il) = iv G Wo div < 2 : v\ E G V(E), VE e E h , Jjvuj dS(x) = 0, Vr G F, J , 

where V(E) = Fq © Pjcc, and |-] r is the jump over T. The density g will be 
approximated in the space of piecewise constants on Eh : 

Q h (Q) = {qe L 2 (n) : q\ E G V£ G . 

Next, we introduce the canonical interpolation operators: 

nf : Wq 1 - 2 n if 2 ' 2 -> : w cur1,2 n w 2 - 2 -> w fc , 

: w div < 2 n w 1 ' 2 - Vfc, n« : l 2 - q„, 

using the available degrees of freedom of the involved spaces. That is, the operators 
(in three dimensions) are defined by [2, 13] 

(nf s) (xi) = sfa), VxieNh] 
J (nfw) x v dS(x) = J^wxi* dS(x), Ve G V h ; 

J (ulv) ■ v ds(x) = J v ■ V ds(x), vr G T h ; 

/ itfq dx= q dx, VE G E h , 

J E J E 
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where Nt it the set of vertices of Eh ■ It is well known that the following diagram 
commutes ([2,3]): 

Wo 1 - 2 n w 2 - 2 Wo cur1 ' 2 n w 2 ' 2 w* w > p n w x - 2 -^-^ Ll 

nf| <j n-j n«j 



Remark 2.7. The interpolation operators 11^, 11^, and 11^, are defined on function 
spaces with enough regularity to ensure that the corresponding degrees of freedom 
are functionals on these spaces. This is reflected in writing W c ' 2 n W 2,2 instead 
of merely W cur1,2 and so on. 

In view of the above commuting diagram, we can define the spaces orthogonal 
to the range of the previous operator, i.e., 

r o, 
h 

rO, 
h 

to obtain decompositions (cf. [2]) 



:= {w h G W h ; curl, w h = 0}^ fl W h , 
V h °' x := {v h G V h ; div x v h = 0}^ n V h , 



W h = V x S h + W%'\ (2.5) 

V h = cuil x W h + v£- L . (2.6) 

The following discrete Poincare inequalities hold [3] 

Kl| i2(n) < C \\div x v h \\ L2{n) , Vv G V^, (2.7) 

IKII^nj^Cllcurl^tBfcll^n), Vw G (2.8) 

where the constant C is independent of h. 

In the subsequent convergence analysis, we make frequent use of the canonical 
projection operators. To bound these we shall need the following ([4,13]) 

Lemma 2.8. There exists a constant C > 0, depending only on the shape regularity 
of Eh and the size of f2, such that for any 1 < p < oo, 



<Ch\\v x 4>\\„ iait 



\\ v - UY lV \\ LP{n} + h\\div x (v - nXv)\\ LP{n) < Ch s \\^ x v\\ LP[n) , r = 1,2, 

\\w - ^fw\\ Lp{n) + h \\cur\ x (w - Il£»|| iP(fi) < Ch s \\W s x w\\ LP{n) , s = l,2, 

for all (j) G W^ p (n),v G W S ' P (Q), and w G W 2 ' P (Q). 

We will also need the following lemma. It follows from scaling arguments and 
the equivalence of finite dimensional norms [ ] . 

Lemma 2.9. There exists a constant C > 0, such that for 1 < q,p < oo, and 
r = 0, 1, 

WhWw.PiE) ^ Ch ~ r T ~ T WMl^E) > 

for any E G Eh and all polynomial functions 4>h G ¥k(E), k = 0, 1, .. The constant 
C depends only on the shape regularity of Eh and polynomial degree k. 

The next result follows from scaling arguments and the trace theorem [1] 

Lemma 2.10. Fix any E G Eh and let <f> G W 1,2 (E) be arbitrary. There exists a 
constant C > 0, depending only on the shape regularity of Eh such that, 

Hh'F) < Ch-? {U\\l HE ) + h\\v x 4>\\ LHE )) , vrer h n dE. 
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We now establish a Sobolev embedding estimate for the discrete decompositions 
(2.6) and (2.5). 

Lemma 2.11. The finite element spaces V?' (CI) and W^'^(f2) satisfies the fol- 
lowing embedding results independent of h: 

(1) The space V^iCl) n W^ v ' 2 (Ct) is embedded in L r {Cl), 

(2) The space w£ ,X (fi) n W curl,2 (f2) is embedded in L 2 * (Cl), 
where 2* = 6 if N — 3, and 2* is any large finite number if N = 2. 

Proof. We first prove (1). By virtue of the decomposition (2.6) we can for any 
v h Gb V,^(n) n W div,2 (O) find functions & € W h (Cl) and js fc G v£' x (fi) such 
that 

li\ (Va-A" 1 [div x Vh]) = curl x £ h + z h . 

Using the commutative diagram and the definition of n^' (V^A -1 [•]) we easily 
verify that 

divz 11^ (V^A -1 [div^ v h ]) = div^ v h . 

Hence, since (zh — Vh) G V?' (Cl) we can use the discrete Poincare inequality (2.7) 
to conclude that 

IK - z h \\ L 2 (n) < C\\ diVaK ^ ^)IIl 2 (o) = 0. 

Thus, Z/j = Vh a.e in f2 and we easily calculate 

IKIIi3*(Q) = \\ z h\\v(n) < ll n ft (V^A^ 1 [divxVh]) \\ L 2* {n) 

< CiUV^A -1 [div^Uh] ||x,a*(n) < C2 1| div K ^11^(0), 

where the last inequality is the standard Sobolev embedding W 1,2 (0) C L 2 (Cl). 

In two spatial dimensions, (2) follows directly from the standard Sobolev em- 
bedding W 1,2 (Cl) C L 2 (Ci). To prove (2) in three spatial dimensions, fix any 
w h G b lf k u (!l)nW„ c " rl ' 2 (!!) and let 77 G W curl,2 (fi) n W div < 2 (ft) C W 1 ' 3 ^) solve 
(cf. [10]) 

curia, 77 = curia, to ft, in fi, 
diva, 77 = 0, in £7, 
77 x f = 0, on dCl. 

Using the decomposition (2.5) of the space Wh(Cl), we can find functions s/j G Sh(Cl) 
and 0» G W^' X (n) such that 

njf 77 = v x s h + Cj.. 

Hence, from the commuting diagram property, we deduce 

curl x Ch = curia; njf 77 = curia- »7 = n ft curia; w h = curia- tOft. 

Thus, since (wh — Ch) £ W®" 1 ^) we can use the Poincare inequality (2.8) to 
conclude that 

\\w h - Ch\\ L 2(n) < C|| curl,,. - Cfc)lk 2 (fi) = 0, 
and hence that Wh — Ch- Moreover, we easily calculate 

IKfelU**(n) = ||Ch||x,3*(n) < I|n^?7lli3*(n) < Ci|Mli**(n) < C 2 \\V x v\\ L^(n) , 

where the last inequality is the standard Sobolev embedding W 1,2 (fi) C L 2 (Cl). 
This concludes the proof. □ 

We end this section by recalling a compactness result from [8, Theorem A.l]. 
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Lemma 2.12. Let {vh}h>o be a sequence of functions in V h °'' L with div x Vh Eb 
L 2 (fl). For any £ G R N , 

\\v h (x) - v h (x - OWmn) < C(\^ + \e)i\\ div x v h \\ L2{n)7 
where the constant C > is independent of both h and £. 

3. Numerical method and main result 

In this section we define the numerical method for the Stokes approximation 
equations and the state the main convergence theorem. The proof of the main 
theorem is deferred to subsequent sections. 

Given a time step At > 0, wc discretize the time interval [0, T] in terms of the 
points t m = mAt, m = 0, . . . , M, where we assume that MAt = T. Regarding the 
spatial discretization, we let {Eh}h be a shape regular family of tetrahedral meshes 
of il, where h is the maximal diameter. It will be a standing assumption that h 
and At are related such that At = ch, for some constant c. Furthermore, for each 
h, let denote the set of faces in Eh- 

For each fixed h > 0, wc let W/j(f2) and Vfc(f2) denote the Nedelec spaces of 
the first order and kind on Eh (cf. Section 2.3) and Qfc(fi) the space of piecewise 
constants on Eh- To incorporate boundary conditions, we let the degrees of freedom 
of Wfe(fi) and Vh(ft) located at the boundary dCl vanish. 

Before defining our numerical scheme, wc shall need to introduce some additional 
notation related to the discontinuous Galcrkin scheme. Concerning the boundary 
BE of an element E, we write /+ for the trace of the function / achieved from within 
the element E and /_ for the trace of / achieved from outside E. Concerning an 
edge r that is shared between two elements E- and E + , we will write /+ for the 
trace of / achieved from within E + and /_ for the trace of / achieved from within 
E-. Here E- and E+ arc defined such that v points from to E + , where v is 
fixed (throughout) as one of the two possible normal components on each edge T 
throughout the discretization. We also write |/] r = /+— /_ for the jump of / across 
the edge T, while forward time-differencing of / is denoted by [/ m ] = f m+1 — f m . 
Discrete implicit time discretization of a function / is denoted by the operator 

Let us now define our numerical scheme. 

Definition 3.1 (Numerical scheme). Let {Qh( x )} h>0 be a sequence (of piecewise 
constant functions) in Qh(£l) that satisfies Q° h > for each fixed h > and g° h — > g° 
a.e. in f2 and in L 1 (fi) as h — > 0. Let the sequence {u°}/j>o be such that for each 
fixed h > 0, u° h G Vh(f2) and satisfies 

/ u° h v h dx= I u°v h dx, Vv h G V h (Q). (3.1) 
Jn Jn 

Now, determine functions 

(Q™> w h n > u 'h) e Qh(ty x Wh(£l) x Vfc(fi), m=l,...,M, 
such that for all <fih G Q/j(0), 



d h t (ft) fa dx = AtJ2 / (Q-« ■ v) + + Q+« ■ ")") I^lr dS(x), (3 2) 
r^r 7 



and for all {-q h ,v h ) G W h (Q) x V h (£l), 

d 1 ? (iC) v h + //curU w h n v h + [(p + A) div,, 11% - p{g^)] &iv x v h dx = 0, 



u 

m M „ , m 

uu 

n 



(3.3) 

w^rjh - Uh curl,, 77^ = 0, 
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for m = 1, . . . , M. 

In (3.2), (tih ■ v) + = maxjit/j • v, 0} and (u^ ■ v) + = mm{ti/, ■ ^, 0}, so that 
uji ■ v = (uh ■ v) + + (uh ■ v)~ , i.e., in the evaluation of g(u ■ v) at the edge F the 
trace of g is taken in the upwind direction. 

Remark 3.2. Recall that g± and (uh ■ v) related to a face T has a different mean- 
ing than g± and (u^ ■ v)^ related to the boundary of an clement dE. By direct 
calculation, one can verify the identity 

At Y, [ (<?+« ' v ) + + Q-W ■ 4>h dS(x) 
E& E h JdE\on 

= - At E / (Q-(<^) + + Q+(<^)~)[<f>h]rdS(x). 
Using this identity, we can state (3.2) on the following form: 




= / Qh'^h dx. 
Jn 

The form (3.4) will be used frequently in the subsequent analysis. 

For each fixed h > 0, the numerical solution {(q™, w™, u™)} m=0 is extended to 
the whole of (0, T) x Q by setting 

(Q h ,w h ,u h )(t) = (Q%,w%,u%), te(t m ^,t m ), m=l,...,M. (3.5) 

In addition, wc set gh(0) = gl and Uh(Q) = u^. 

The continuity scheme (3.2) clearly preserves the total mass. The following 
lemma from [8, Lemma 4.1] states that the density is strictly positive whenever the 
initial density is strictly positive. 

Lemma 3.3. Fix any m = 1,...,M and suppose g™^ 1 G Qh(Q), u™ G Vh(fl) 
are given bounded functions. Then the solution g™ G Qh(ty of the discontinuous 
Galerkin scheme (3.2) satisfies 

minp™ > ming" 1 - 1 ( „ - 1 n ) . 

Consequently, if g™^ 1 ^) > 0, then g™(-) > 0. 

Existence of a solution to the nonlinear-implicit discrete scheme follows from a 
topological degree argument. This argument is essentially identical to that of [8, 
Lemma 4.2] with a minor modification to accommodate the discrete time derivative 

Lemma 3.4. For each fixed h > 0, there exists a solution 

(Qf,w^,uf)eQ h (Q)xW h (Q)xV h (n), £/,"(•) >0, m = l,...,M, 

to the nonlinear -implicit discrete problem posed in Definition 3.1. 

Our main result is that, passing if necessary to a subsequence, {(gh,Wh, Uh)} h>0 
converges to a weak solution. More precisely, there holds 

Theorem 3.5 (Convergence). Suppose (g°,u°) G n I 2 (SJ), 7 > f ■ Let 

{(gh,Wh,Uh)} h>Q be a sequence of numerical solutions constructed according to 
(3.5) and Definition 3.1. Then, passing if necessary to a subsequence as h — ► 0, 
Uh — ► u, a.e in (0,T) x Q, ghUh — gu in the sense of distributions on (0,T) x Q, 
and gh — ► g a.e. in (0,T) x Q, where the limit triplet (g,w,u) satisfies the mixed 
formulation (2.3), and thus (g,u) is a weak solution according to Definition 2-4- 
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4. Basic estimates 

In this section we gather some basic estimates for our numerical method. The 
results include stability and weak time-continuity of both the density and velocity. 
We however commence by recalling (from [8] ) the following renormalized version of 
the continuity scheme. 

Lemma 4.1 (Renormalized continuity scheme). Fix any m = 1, . . . , M and let the 

pair (g™,u™) G Qh x Vh satisfy the continuity scheme (3.2). Then (q™,u™) also 
satisfies the renormalized continuity scheme 

B(Qh)<t>h dx 

-MY, I (B(Q m )« ■ v) + + S(^ l )« • f)") Wr dx 

+ At [ b( Q %)div x u^ h dx+ I B"{H{Qk,Qh- 1 )) bh' 1 } 2 ^ dx 

Jn Jn (4.1) 

+ At Y I B"(a T {e7 . Q m )) iQhfv ■ ^) + 

- B"(?(g™, ($)) [e^Jr (M+(uf ■ v)- dS(x) 

B(g™~ 1 )4>h dx, V0 fc G Q h (n), 

for any B G C[0,oo) n C 2 (0,oo) with B(0) = and b{g) := qB'(q) - B(g). Given 
two positive real numbers a\ and 02, we denote by £((11,02) and £ r (eti,a2) two 
numbers between a\ and ai (See [8] for a precise definition). 

In what follows we will need the following discrete Hodge decomposition. 

Lemma 4.2. Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. For each fixed h > 0, there exist 
unique functions £™ G Wi' and z™ G V^ ' such that 

u% = curl x Ch+ z h, m = 0,...,M. (4.2) 
Moreover, if we let Ch(t,x), Zh(t,x) denote the functions obtained by extending, as 
in (3.5), {C}m=i> {*r>m=i to </le whole °/( ; T ] x </le71 

u h (t,-) = eurl x Ch(;t) + z h (;t), te (0,T). 

Finally, let curl, £° G Z 2 (f2) and V,s° G i 2 (r2) satisfy the standard continuous 
Hodge decomposition u° = curl, £° + V,s°. Then, 

curl, C° h -> curl, C°, z£ -> V K s°, in L 2 (ft), 

where Ch an d are given by (4.2). 

Proof. The first two statements are consequences of (2.6). 

To prove the last statement, fix any <fi G C£°(fi) and set Vh = 4> m (3-1) to 
obtain 

/ curl, C° h curl a (nf 0) dx = / curl, C° curl,(nf <p) dx, (4.3) 

where we have used that it° = curl, £° + z\ and / n z° curl, II dx = 0. Now, 
since || curl, Chilian) < C|l u ° lli 2 (o) < C|| w° || i2(n) , there exists a function curl, C° 
such that curl, £° — ^ curl, £° in L 2 ( fl). Sending h — * in (4.3) yields 

/ (curl, C° - curl, C°) curl, dx = 0, V0 G C c °°(fi). 
Jo 
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Hence, curl, £° = curl, £° a.e in O. 

Next, let Vh = curl, in (3.1) to discover 



curl^ Chilian) = J curl * £° curl * Ch dx -* II curl * C 



|| 2 

£ 2 (n)' 



as h 0. Then, curl, £° -> curl, C° in L 2 (f2). 
By setting = u° in (3.1) we deduce 



l«°lll=(Q)= / uM<fc-HI« lll»cn), 



as ft — ► 0. Hence, — ► u° in L 2 (f2). 
Finally, a direct calculation shows that 



„0||2 



= lim ||u° - «°||i a(n) - lim || curl, $ - curl, C°lli»(n) + \\4 ~ V * s ' ll£»(n 



- 2 lim 

h.->0 



/ (curl, Ch - curlx C°) (*° - V,s°) 



where the last term converges to zero since curl, £° — > curl, £° in L 2 (f2). Thus, 
z° — > V,s° in L 2 (f2) and the proof is complete. 

□ 

We now derive a basic stability estimate satisfied by the numerical scheme. 

Lemma 4.3 (Stability). Let {(g n , Wh, u h)}f l> o be a sequence of numerical solutions 
constructed according to (3.5) and Definition 3.1. For g(-) > 0, let 

7 — i z 

For any m = 1, . . . , M , there holds 

/ E{el\K) dx + Y, At|mjj 2 wdiv , 2(0) M\\w k h \\ 2 w ^, 2{Q) +K 

< f £(g°,u°) dx, 
where the numerical diffusion term is given by 



diffusion 

fc=l fc=l 



J v diffusion ^ 



I I L lib r. 

K _1 1 HW>)+£ / P "^- h (QlQ k h - 1 )){Q k h - l f dx 

k=l k=l Jn 
m „ 

+ EE At P"{g k ,){g k h f T K-v\ dx. 
fc=irer£ r 

Proof. The proof is almost identical to that of Lemma 5.3 in [8] and follows directly 
from standard arguments. We omit the details. □ 

Since the finite clement spaces are not conforming in W 1,2 (£l) it is not clear that 
the velocity and vorticity are embedded in L? (f2). Knowing this is essential for 
the later convergence analysis. 

Lemma 4.4. Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Then 

w h e 6 L 2 (0,T;L 2 *(r»)), u h e b L 2 (0,T; L 2 *(f!)), 

where 2* = 6 if N = 3 and 2* is any large finite number if N = 2. 
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Proof. The second equation in (3.3) with test function rjh = V x s/j reads: 

/ w™V x s h dx = 0, Vs h £ S h (Q), m = 1, . . . , M, 
Jo 

where the space Sh(Ct) is defined in Section 2.3. By definition, this means that wV? £ 
Wl' (0) and hence Lemmma 2.11 is applicable and yields the desired estimate: 

w h £ h L 2 (0,T:L 2 *(tt)). 

Next, we make use of Lemma 4.2 and let {Ch}h>o, {zh}h>o satisfy 

u h (-,t) = cm-l x Ch(-,t) + Zh(-,t), 

C h (-,t)£W°' ± (ri), z h (;t)ev°' ± (a.), 

for all t £ (0,T). 

Another application of Lemma 2.11 yields 



(4.4) 



z h £ h L 2 (0,T;L T (Q.)). 
Fix »7 e Wq^'^^^O) and let rj h £ W°^(ft) satisfy 



(4.5) 



/ curl x rjh curia <f> h dxdt = / curl x rj cur\ x dt h dxdt, \/<j> h £ W h (Cl). 
Jn Jn 

Then, by utilizing the second equation in (3.3) with rjh as test function (the second 
equality below), we calculate 



M 



curL rj cur\ x £™ dx 



M . 

E A * / < 



rffdxdt 



M . 

V] At / curia; r/,, curls ^ 

m=l 



< H' u 'A||jr2(o,TjIi a *(n))ll T W>ll£(a*)'(n)) 

< II' u, 'i!Il 2 (o,T;I, 2 *(0))II cur ^ , 7lli(2*)'( n ) I 



(4.6) 



where the last inequality follows from the discrete Poincare inequality (2.8). 

Now, for an arbitrary d> £ hS 2 > let curL r\ be given through the Hodge decom- 
position <f) = curLj rj + X7 X \. Then, we can use (4.6) to deduce 

\f Q d>cur\ x C h dx\ 



sup 11^,11 



dt 



sup 

. 0en 2 *)'(n) 



sup 

o \0eL< 2 *>'(fi) 



| J n cur^ r} curia, Oidx| 



I0lli( 2 *)'(o) 



\Sq VhW h dx\ 
ll<AII.L( 2 *>'(n) 



dt < \\wh\ 



i 2 (0,T;I, 2 *(n))' 



where the last term is bounded from (4.4). Hence, curb; Ch £b L 2 (0,T; L 2 (CI)) 
and, keeping in mind (4.5), u h £ b L 2 (0,T; L 2 * (0)). □ 

In the upcoming convergence analysis and in order to establish weak time- 
continuity of the density we shall need to control the artificial diffusion introduced 
by the upwind discretization of the continuity equation. The following lemma pro- 
vides the required bound. 

Lemma 4.5. Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Then there exists a constant C > 
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depending only on the initial energy E (g , u ) , the shape-regularity of Eh, T, and 
|0|, such that 

(■T 



E 

EeE h 



IqhI (u h -v)-(n%<t>-<j>) dS{x)dt 



'0 JdE 

</i e(7) C||V^|L 2(0>x . L2 . (fi))! V0eL 2 (O,T;W^ 2 >)), 

where 2* = 6, if N = 3, and 2* is a sufficiently large number, if N = 2. Here, 
0(j) > is given by (4.14) below. 

Proof. Let 4> e £ 2 (0,T; W 1 ' 2 *^)) be arbitrary and set 



1 

At 



(s,x)ds, C = n^ m , m = l,...,M. 



We will need the auxilary function B(z) = z a . where 



7 



and i £ N is chosen such that 76 (i + 1, i + 2]. 



i + 1 

Using B"(z) > and the Holder inequality, we obtain 

A/ . 



< 



E E A * 

m=lE£E h JdEXdQ. 
/ M 

E E At 

\m=lE£E h JdE\dfl 
( M 

x E E^ 



B"{ Q ™)iQ™f\u™-v\dS{x)j 



h x I 2 - 



dE\dn 



In the case ^<7<2, a = 7 and Lemma 4.3 yields 



h<C [ B(qo) dx = C [ (g°y dx. 
Jn Jn 



Conversely, if 7 > 2 then 2a < 7 and the renormalized scheme (4.1) with <fih := 1 
yields 

M 



h<{a- 1) 



E Ai / (eD^div.Mr^ + / (f?T^ 



< C (llM2«>(0,T;Z/?(O))ll div z «ft|U 2 (0,T;£2(fi)) + J^Vdx 

which is bounded by Lemma 4.3. Consequently, in both cases, we conclude that 

h < C. (4.7) 
To bound the I2 term, we utilize the Holder inequality: 



\fgh 2 \ dxdt 



Jn 



< 



\fg\ 2 *- 2 dx 



\Jn 

T 



\hf dx I dt 



(4.8) 



< 



I/I mi dx) \\g\\ L 2- ( n)\\h\\b- { n) dt 

\JQ / 

< ll/!|L^(0,T;L m i(fi))|Mil,=°(0,T;,L 2 *(n))IH! L 2 (0,T;L 2 * (fi)) ' 
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where I <m\ = tt|^t < 2 and — + i + £ = 1. 
Now, using (4.8), we deduce 



\EeE h 



M 



*EM E 

m=l \EEE h 



dE\dn 



n; 



2' 



J 2 <a(o-l) max [E / ' M ™ ' "I* d5 ( a; ) ) 

V 

dS{x) 

(B"(g?)) 1 dS(x) 



(4.9) 



x max 

m=l,...,M 



E 

\E£E h 



!dE\dCl 

Next, we apply Lemma 2.10 to deduce 



max 

m=l,...,M 



E 

\EdE h 



ldE\dn 

Similarly, we find that 

M 



K-H 2 * dS{x)\ <Ch-^\\u h \\ L ^ T . L .» m . (4.10) 



E a m E 

m=l \E£E h 



dE\dQ 



n. 



2 ! 



<ch ^||n^-0|| 2 L2(OT;L2 , (n)) 



dS(x) 



< Ch 2 



^llL 2 (0,T;L 2 *(O))' 



where the last inequality is an application of Lemma 2.8. 

To derive a similar bound for the B" term in (4.9), we first note that, since g™ 
is everywhere positive and 2 — a < 1 , 



{b"(q?)) 1 



< 



(2— a)mi 



<C(l + |£™|" ll + |e m | 



on every F n <9i? \ <90. From this, we conclude that 

dS(x) < h- x c \E\ + 



8E 



M{E)ISE 



\Qh\ mi dx , 



where Af(E) denotes the union of the neighboring elements of E. Applying this 
together with Lemma 2.10, we obtain 



max 

,=1,...,M 



_ — _ P i mi 

E / (B"(QT)) dS(x) 

^c. JdE\dCl 



\EeE h 



<ch ™i nn| m i + ||£»/iIIl-(o,t ; l™i(o)) 

, / m i n {0,A r (^--i)} ii || 



(4.11) 



where the last inequality is a standard inverse estimate (Lemma 2.9). Setting 
(4.10)-(4.11) into (4.9) leads to the bound 

h < C^-!l w /i|lL~(O,T;I, 2 *(O))ll^ 7 K0|li2( OiT;L 2*( n ) ) 

(i , ; min{0,Af(— --}\\ || \ 
X(l + ft Vmi T r ||^|U«>(0,T;i-r(n)) J 

i 2 (4-12) 

< C/l 2 ||Mh||L 2 (o,T;Z, 2 *(O))l|V a: 0|| L2( ,T;L 2 *(H)) 

/-, , , min{0.JV( — --}n II A 
X[l + h W "i ^ S \\Qh\\L^(0,T;L-r(n)) j , 
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where the last inequality is an application Lemma 2.9 in time (keeping in mind 
At = nh). We have also used that 

In 2D, 2* is any large finite number. Consequently, we can always make sure 
that mi < 7. Using this, it is straight forward to check that 

where 

f i N = 2 

O<0( 7 ):= ln „ (4.14) 

+ mm{0,3(± - iV = 3, 

By setting (4.13) into (4.12) and applying Lemma 4.3, we obtain 

h<h^C\\VMU ,T,L^ m - 

This and (4.7) gives 

I 2 = h x I 2 < C\\W x cf>\\LHo^ m h 2eh \ 
which brings the proof to an end. □ 

4.1. Weak time- continuity estimates. We end this section by establishing weak 
time-continuity of the density and velocity. 

Lemma 4.6. Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Then 

d?(g h )e b L 2 (0,T;W-^ 2 *y(n)), 

where (2*)' = 2 «_i an d 2* is as in the previous lemma. 

Proof. The proof is almost identical to the proof of Lemma 5.6 in [8] and is only 
included for the sake of completeness. 

Fix <f> £ L 2 (0, T; W 1 ' 2 (f2)), and introduce the piecewise constant approximations 

<f> h := V$<f>, 4% := IL^"\ and 0™ := ^ <f>{t, ■) dt. 

The continuity scheme (3.2) with <\>™ as test function reads 

At / d?(g%)<t> m dxdt 
Jn 

^ r (4-15) 
= A* £ / ■ ^) + + Q+« ■ vT) IClr dS(x). 

rerf/ r 
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Since the traces of 4> m taken from either side of a face are equal, we can write 
E / (^«^) + + ^«^)-)[Clr dx 

• u) + + g™(K ■ v)~) (C - <T) 

E / -div x (^<(c-0 m ))^ 



rer; 



E 

EeEh JaE\aa 



E£E h - 



E 



dE\di1 



in EeEh JdE\an 



(4.16) 



To conclude the last equality, we have used 



e£ div, <(C - 4> m ) dx = ( e £ div, / n^ m - 0™ = o, V£ e E h , 

J E 

since both and div^ it™ are pieccwise constant. 

By summing (4.15) over m, taking absolute values, and using the above identity, 
we find 

M 



V At / dHQh)4> m dxdt 

M 

E At / qZKv,,*" 1 dx 

m=l 



E E A * 



m=l E£Eh JdE\dCl 

Using Lemma 4.5, together with an application of Holder's inequality, we deduce 



m . 

E A< / d' t l (gT)<P m dx 



m— 1 



lL 2 (0,T;i 2,, (n)) 



< E AtlleriU-wlKILa-cnjIIVx^lli^^ + ^wilv^ 

m— 1 

< II £/i || L°°(0,T;L"(fi)) II Wh||x,oo(o,r;i3*(n)) II ^</>|| L 2 (0,T;I, 2 * (Q)) 

+ c^||v^|| L2(0:T;i2 . (n)) , 

where a = 2 «-(2»)' ^ 7 smcc 7 > T an d „ + + ^ = 1- By Lemma 4.3, the 
right-hand side is bounded, so we conclude that 

dt {Qh)^ dxdt 



At Jn 

M 



E At / dt( e z)<r dx 

m=l ^ 



<C(1 + ^^)||V^| 



L 2 (0,T;i 2 *(n)) 
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Jn 

T 

h 



This brings the proof to an end. □ 

Lemma 4.7. Let {(gh, Wh, u h)}h>o ^ e a se Q uen ce of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Then 

d^{u h )& b L 2 (f),T;W- 1 ' l {n)), 
Proof. By adding and subtracting, we see that for any <fi € L 2 (0, T; Wg 1 ' 00 ^)), 

d£ (uh) <f> dxdt 

d% (u h ) U\<t> dxdt + [ [ (u h ) (0 - Iljf 0) dxdt. 
o Jn Jo Jn 

From the first equation of the momentum scheme (3.3) with Vh = nYtfi, we have 
that 

/ / df (u h ) TlX&dxdt = - I n cur\ x w h (njf 0) + (jj, + A) div x u h A\v x H%4> dxdt 
Jo Jn Jo Jn 

+ / a q1 div x nj^0 dxdt 
Jo Jn 

< C (|| CUr^ W/j£2(0,T;£2(n)IMU2(0,T;£2(n) 

+ || diVa; Mft||L 3 (0 : T;L 2 (f2))ll div^ 0| \ L 2 ( ,T;L 2 (Q)) 

+ \\Qh\\L^(o.T-L-<(n))\\ div x 0||ii(o,r ; i~(n))) 

< Cll^llia^.TjW 1 - 00 ^))! 

where the last inequality follows from Lemma 4.3. 
From Lemma 4.3, we also have the estimate 



\\dHu h )\\L H o,T;L H n)) = m-i\^JjK' 1 f dxj <h~?C. (4.17) 
Using (4.17), we estimate 

[ T [ d(n T fe) {4> - nl0) dxdt < c\\d? (u h ) \\mo,T; LH n))U - ul4>\\ LH0 .T;mn) 

J At Jn at 



< C-^=\\V x (t>\\L^(0,T;L2(n)) < Ch^ 

v At 



where we in the last inequality have used the relation At = nh. Combining the 
previous estimates concludes the proof. □ 

Recall our notation for the Hodge decomposition of the solution u, 

u = curLj £ + \7 x s. 

In the next result, we prove that dt(cm\ x £) E L 2 (0, T; L 2 (fl)). To see why such a 
bound is reasonable, apply the curLj operator to the velocity equation (1.5) 

curL(curL £) t + [i curL curL w = 0. 

Multiplying with £ t , integrating by parts in space, and applying Holder's inequality, 

C 

II cur^ Ct\\h(Q) < 4 curia, Ct|li2(n) + — II curl * w \\h(n)- 
Fixing e small, and integrating in time 



/ II cur1 ^ Ctllivn) dt<C [ || curia, w\\ 2 L2(n) dt, 
Jo Jo 
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where the right-hand side is bounded. Consequently, it is the higher regularity on 
w = curlr u that enable us to obtain the bound. 

Lemma 4.8. Let {(gh, w^, Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Let {(C/u z h)}h>o ^ e the sequence 
given by the decomposition Uh(-, t) = curl x Gh(',t)+Zh{',t) and Ch(',t) <= WP' (f2), 
z h (-,t) E V r h °'- L (n), forte (0,T). Then 

$ (curl, Ch) e b L 2 (Q,T;L 2 (n)). 

Proof. For any m = 1, . . . , M, let = 9j (curl x £™) 6 V/j. Observe that by the 
orthogonality of the Hodge decomposition, 

^«)^(curl x CD dx= [ (^(curl.CDl 2 dx. 

Hence, by setting v™ as test function in the first equation of the momentum scheme 
(3.3), multiplying with At, and summing over all m = 1, . . . , M, we obtain 

M . 

VAt / | (cm-LCD | 2 dxdt 

m=l Ja 

M 

= - V At / ncmlxW^d? (curl^CD dxdt 
m=l ^ 

(Af \ 5 / M \ 2 

£ At||d t ' 1 (curl^n ii 2 E^ii^^nii^n) ■ 
m=l / \m=l / 

An application of the Cauchy inequality with e to the above estimate yields 

\\8t (curia, Ch) ||i2(0,T;£a(n)) < ^ II CUrLj Wfc || L2( 0)T . i 2 (n )) . 

Lemma 4.3 provides a bound on the right-hand side and hence the proof is complete. 

□ 



5. Higher intergrability on the density 

The stability estimate only provides the bound p(gh) Gb L°°(0, T; L 1 (fi)). Hence, 
it is not clear that p{gh) converges weakly to an integrable function. Moreover, the 
subsequent analysis relics heavily on the pressure having higher integrability. In 
this section we establish that the density is in fact bounded in L 7+1 (0, T; L 7+1 (f3)), 
independently of h. The main technical tool used to achieve this is an equation for 
the effective viscous flux: 

P e s(g, it) = p(g) - (A + (i) div x u. 

We start by deriving this equation. For this purpose, fix any <f> <G L 2 (0, T; Lg(f2)) 
and, for each fixed h > 0, let 

v h (t, ■) = III (v x A- 1 [</>]) (*,•), te(0,T), 

and 

1 /•*" 

< = XtJ , V/l(s ' ^ ds ' TO = : ' ■ ■ ■ ' M - 

Observe that is constructed such that 

1 /■*" 

div x v% = y i dfc, m = 1, . . . , M. 
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By inserting v™ as test function in the momentum scheme (3.3), multiplying 
with At, and summing over all m = 1, . . . , M, we are led to the identity 

f / P e s(Qh,U h )(/) dxdt= [ [ (d^( Uh )+fJ,CUTl x W h )UX (V^A" 1 [0]) dxdt. 

o Jo Jo Jo 

Since Jo(curLj w™)V x A -1 [<f>] dx = 0, for all m = 1, . . . , M, we can further write 



P c s{Qh,u h )4> dxdt 



o Jo 

T 



J!) 



/ ^ (tt h ) V^A" 1 [0] dxdt 
Jo 

+ [ f (u fc ) + /i curl, «7 h ) (itf (V, A" 1 [</»]) - V, A" 1 [<£]) dxdt. 
Jo Jo 

(5.1) 

As ^ was fixed arbitrary, we can conclude that (5.1) holds for all <f> £ L 2 (0, T; Lq(SY)). 
The following lemma ensures that the last term of (5.1) converges to zero. 

Lemma 5.1. Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Then there exists a constant C > 0, 
depending only on the initial data and the shape regularity of Eh, such that 



f [ (u h ) + M curL, w h ) (l# (V^A- 1 [<f>]) - V^A" 1 [ft) dxdt 
Jo Jo 

< C(hi + h)M\\ L '<f>,T;L'<n))> e l2 (°' T ' L o(^))- 
Proof. By this, the Holder inequality, and Lemma 2.8, we deduce 

/ / (d? {u h ) +fxcuvl x w h ) {TLX (V.A- 1 [0]) - V.A" 1 [0]) dxdt 
Jo Jo 

< c/lUV^V^A^ 1 [<f>] \\l*(0,T;L»(C1)) 

X (||a t ' 1 (Uh) \\l'(0,T;L'(SI)) + II CUrL B 1i>/,||£2(o ) T;£a(n))) 

where we in the last inequality have used Lemma 4.3 and (4.17). 



□ 



We are now in a position to prove higher intcgrability of the density. To increase 
readability of the proof, we introduce the notation 

for the spatial average value of a function. 

Lemma 5.2 (Higher intcgrability on the density). Let {(gh,Wh,Uh)} h>Q be a se- 
quence of numerical solutions constructed according to (3.5) and Definition 3.1. 
Then 



Qh e 6 LT +1 ((0,T) x SI). 
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Proof. Setting <f> = £>/, — (gya in (5.1) yields the identity 



np{Qh)Qh dxdt 
A 



JQ 



p{Qh)(Q° h )n + {>> + p)div x u h g h + d? {u h )V x A 1 [(g h - (qDq)] dxdt 



+ / (pi (u h ) +fj,curl x w h ) 
Jo Jn 

x (TLX (V.A- 1 [g h - (g° h ) n ]) - V^A" 1 [g h - (g° h ) Q ]) dxdt, 
Applying the Holder inequality and Lemmas 4.3 and 5.1 yields 



p(gh)gh dxdt 



o Jn 



< 



I f flf (waJVxA- 1 [g h -(g° h )n] dxdt 
Jo Jn 

+ C (l + + h) \\gh\\Li(o,T;L2(n))- 



To bound the first term on the right- hand side, we first note that 



(5.2) 



o Jn 



d£ (u h ) V^A" 1 [g h -(g° h )n] dxdt 



J2 At / V.A- 1 [<ff - <<£>„] dx. 



(5.3) 



Then, we apply summation by parts to (5.3) and make use of the Holder inequality 
to obtain 



(5.4) 



/ / d£ (u h ) V X A _1 [g h - (g° h )n] dxdt 
Jo Jn 

M r 

-J2 At uJT'VsA- 1 (#)] dx 
k=i Jn 

A/ . 

]TAi / uJ^VaA" 1 [#(#)] dx 

m=l Jn 

+ C|l u °lli, 2 (n)ll£ l /i - (gh)n\\L°°(o,T;L-y(n)), 

where we in the last inequality have used elliptic theory (and 7* > 2 since 7 > -y) 
to conclude that 



\\7 X A 1 [g h - (g h )n] ILoo( 0i T;£2(n)) - c \\Sh ~ (fti)n|U°°(o,T;£T(n))- 
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Next, using integration by parts, 

£>t / uJT^A" 1 dx 

, Jo. 



m—1 

M 



J2^t / A- 1 [div^/- 1 ] 



At 



dx 



V At f m^^V.A- 1 [div.M™" 1 ] dx 

m=l Jn 



M 



+ E At E / IAe « ■ - I) A" 1 [div, u™- 1 ] dS(x), 

m =l E^E h JdE 

(5.5) 

where the last equality is deduced as follows: Set <\>™ = 11^ A -1 [div^ in the 

continuity scheme (3.2) and perform the calculation (4.16). 

By setting (5.5) into (5.4), and applying Lemma 4.5, we obtain 



Jo Jq 



Qh 



{Qh)n\ 



dxdt 



< C 1 



|ttfc|l L»(0 I T; i ^ r (n)) llefc||iOO(0 ' Tii ' ,r{n)) 



(5.6) 



+ /i fl WC||V 1B A- 1 [div x u h ] || L2(0 , T;i2 . (o)) , 

where 0(7) is given by (4.14). 

Finally, inserting (5.6) into (5.2) and recalling that < 2*, since 7 > -y, gives 

T f 1 
Jn 

< C(l+ ||Mh||L 2 (0,T;i 2 *(O))lle'i||L° c (0,T;LT(n))) + ^ 7) C|| div^ Mh||L 2 (0,T;L 2 (n)) 
+ (l + ^ + ft) ||£ l /l||L 2 (0,T;L 2 (O))- 

The proof then follows from the Holder and Cauchy (with epsilon) inequalities. □ 

6. Convergence 

Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions constructed accord- 
ing to (3.5) and Definition 3.1. In this section we establish that a subsequence 
of {(Qh, Wh, Uh)}h>o converges to a weak solution of the semi-stationary Stokes 
system, thereby proving Theorem 3.5. The proof is divided into several steps: 

(1) Strong convergence of the velocity 

(2) Convergence of the continuity scheme. 

(3) Weak sequential continuity of the discrete viscous flux. 

(4) Strong convergence of the density. 

(5) Convergence of the velocity scheme. 

Our starting point is that the results of the previous sections ensure us that the 
approximate solutions (wh,Uh, Qh) satisfy the following ^.-independent bounds: 

Qh G b i oo (0,T; J L''(r!))nLT+ 1 ((0,T) x fi), 

w h e b L°°(0, T; L\n)) n L 2 (0, T; W curl ' 2 (ft)), 

u h e b L°°(0, T; L 2 (ty) n L 2 (0, T; W^W). 
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Moreover, in view of Lemma 4.2, there exists sequences {Ch}h>o, {zh}h>o such that 
u h (-,t) = curl, Ch(-,t) + z h (-,t), 

a(-,t)e< A (^), z h (- t t)eV°' x {Sl), 

for all t e (0, T) where 

z h e b L°°(0,T;i 2 (r!))nL 2 (0,T;W div ' 2 (fi)). 
curl, Ch £ h L°°(0,T;L 2 {Q)), 

and 

# (curl, Ch) Gb L 2 (0,T;L 2 (Q)). 
Consequently, we may assume that there exist functions g,w,u such that 

0h g, in L°°(0,T;LT(fi)) n^((0,T) x 0), 

tflfc^to, inL oo (0,T;i 2 (fi))nL 2 (0,T;W curl ' 2 (fi)), (6.2) 

« fc it, in L°°(0, T; £ 2 (ft)) n L 2 (0, T; W div ' 2 (O)). 

Furthermore, using the standard Hodge decomposition u = curl, C + V,s and 
orthogonality, 

z h -°V x s, inL oo (0,T;i 2 (fi))nL 2 (0,T;W div ' 2 (fi)), 

(6.3) 

curl, Ch curl, C, in C(0, T; Z 2 (ft)) n V^ 1 ' 2 (0, T; £ 2 (ft)). 
In addition, 

0ft P 7 , 0ft P 7+ , 6h log Qh QlogQ, 

where each h ^-° signifies weak convergence in a suitable LP space with p > 1. 

Finally, p h , g h \ogg h converge respectively to g, q log gin C([0, T]; L^ eak ( »)) for 
some 1 < p < 7, cf. Lemma 2.2 and also [5,11]. In particular, g, glogg, and g\ogg 
belong to C7([0,r];i^ eaj£ (O)). 

6.1. Strong convergence of the velocity. 

Lemma 6.1. Let {(gh,Wh,Uh)} h>0 be a sequence of numerical solutions con- 
structed according to (3.5) and Definition 3.1. Then 

u h ^u, in L 2 (0,T;L 2 {Q)). 

Proof. By virtue of (6.1) we can consider each component of the decomposition 
Uh = curl, C/i + Zh- In Lemma 6.2 below we prove that 

curl, Ch curl, £ in L 2 (0,T: i 2 (fi)), 

and hence it only remains to prove that Zh z m the sense of distributions. 
From Lemma 4.7, we have the the weak time-continuity estimate: 

d*{z h )£ h L 2 {Q,T-W- l >\n)). 

Lemma 2.12 provides the spatial translation estimate: 

Q 4— N 1 

\\z h (t,x) - z h (t,x- OIU 2 (o,T;£ 2 (n)) < C(ICI + l£l~) 2 ll div x z/ l || i 2(o i T;i,2(a)) J 

where the constant C > is independent of h and £. Lemma 2.3 can then be 
applied (recalling (6.3)) to obtain the desired result; 

z h ^V x s, mL 2 (0,T;L 2 (n)). 

□ 
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Lemma 6.2. Given (6.2) and (6.3), 



w h w, curl, Ch curl, C in L 2 {0, T; L 2 {Q)). 

Proof. Fix any * G (0, T) and m such that t G {t m - x ,t m \, where t m = mAt. 
Subtract the first equation of (3.3) with Vh = curl, £™ from /x times the second 
equation of (3.3). Multiplying the result with At and summing over all k = 1, . . . , m 
yields 



h->0 



o Jn 



H curl, r] h curl, Ch - /z curl, curl, £/, dxdt 

: / / izw^rj/j + <9 t ' 1 (curl, Ch) curl, & dxeft, 
Jo Jn 



(6.4) 



for all r]h,^h that are piecewise constant in time with values in W/j(f2). Fixing 
77, £ G C£°((0,T) x Q), we use in (6.4) the test functions 

1 /•*'" 

rih{t,-)=Vh (■)■■=-£. J m K h v n {;s)ds, t€(t m - 1 ,t m ),m = l,...,M. 
^(*'-)=C('):=^/* ^^(..sjds, te(Vi,g,m = l,...,M. 

Due to Lemma 2.8, curl, £/, — ► curl, £ and curl, rj^ — > curl, 77 in L 2 (0,T; L 2 (Q)). 
As a consequence, keeping in mind (6.2) and (6.3), we let h — ► in (6.4) to obtain 



Jn 



fx curl, 77 curl, C — M curl, u? curl, £ cfedi 

: / / fiwr) + d t (curl, C) curl, £ dxdt, V77, £ G C c °°((0, T) x fl). 
Jo Jn 



(6.5) 



Since C c °°((0,r) x SI) is dense in i 2 (0, T; W curl *' 2 (ft)) ([10]), we see that (6.5) 



holds for all 77, £ G £ 2 (0,T; W curll ' 2 (fi)). Hence, taking 77 = u>, £ = C in (6.5), 

curl,C°| 2 dx = / / /x |r«| 2 dzeft + ^ ( /" |curl,£| 2 dx] (t), (6.6) 



where curl, C° is given by the Hodge decomposition it = curl, C° + V,s°. 
Next, setting rjh = and £^ = Ch in (6-4), we observe that 

- / I curl, C"| 2 dx = I [ fJ.\w h \ 2 dxdt 
z Jn Jo Jn 

+ \ ( 1 1 cmi, a i 2 dx) (t) + \ Yl [ i curl - ^"i 2 dx - 

\Jn J m _ x Jn 



(6.7) 



Subtracting (6.6) from (6.7) 

/ / fi(\w h f - H 2 ) dxdt curl, C/,1 2 - I curl, C| 2 dx) (t) 

Jo Jn * \Jn / 

±J I curl, C£| 2 - I curl.C ! 2 dx 



lim 



< lim 
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where the last equality is an application of Lemma 4.2. Consequently, for any 

*e(o,T), 



lim 



fi\\Wh - W\\ L 2 

t 



,t;L*(Q)) + \\\ curl * Ch(t, •) - curls CO, Olll^n) 



= lim 



/ (i (\w\ 2 — Whw) dxdt 
Jo Jn 

+ Dm (J | cur^ C| 2 - (curl K GXcurL, C) dxj (t) = 0, 
where the last term converges to zero due to the weak convergences (6.3). 



□ 



In the subsequent analysis, we will need the following technical lemma. For 
notational convenience, we define the linear time intcrpolant Tic: 



(U c f) (t) = f m - L + 



At 



(f m -f m - L ), t€(t m -\t m ). 



Lemma 6.3. Given (6.2) and (6.3), 

A- 1 [div x u h ] A-^diVsu], in L 2 (0, T; W 1:2 (fl)), 
A-^div^cUh] A-^div^n], m tfiO.T-W 1 ' 2 ^)). 
Proof. Recall the continuous Hodge decomposition 

u = curlj; £ + VsS. 

As in the proof of Lemma 2.11, we have that Zh = (V K A _1 [div^ Uh]). Hence, 
||n^ (V^A" 1 [div x Uh\) - V X A _1 [div x u] \\l 2 (o,t-l 2 {Q.)) 

= \\Zh - V a! s|| i 2(o,T;I, 2 (f2))- 

From Lemma 6.1, we have that the right-hand side converges to zero. Hence, we 
conclude that li\ (V^A" 1 [div x u h ]) -> V^A^ 1 [div^ it] in L 2 (0, T; L 2 (Vt)). Next, 
we write 

||V 2 ,A _1 [div x Uh] - V^A" 1 [div x u] \\ L 2 (o.T-L 2 (n)) 

< HV^A -1 [div x u h ] -11% (V^A -1 [div x u h ]) \\L*(p,T;L'(n)) 

+ \\HX (V^A- 1 [div x u h ]) - V^A" 1 [div x u] \\ L *(o,T;ma)) 

< Ch\\ div x Mh||L 2 (o,T;L 2 (n)) 

+ (V X A^ [div x u h ]) - V^A" 1 [div x u] Wvqmww 

By sending h — > 0, we discover 

V.A" 1 [div, u h ] -► V.A" 1 [divx u] in L 2 (0, T; L 2 (fl)), 

which proves the first part of the lemma. 
A direct calculation gives 

\n zUh {t r )-u h {t,-)\ a <\{u^H-)} 2 \, te(t k -\t k ). 

Hence, integrating over (0, T) x Q yields 



M—l 



\\Iicu h - u h f L2(0tT . L2m < At || Kl || i2( n) ^ CAt > 



(6.9) 



k=l 



where the last inequality follows from Lemma 4.3. 
Using elliptic theory and (6.9), we conclude 



I A" 1 [n c u h -u h ]\ 



L 2 (0,T;W 1 . 2 (n)) 



< C(Ai)3. 
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Hence, the limits are equal and consequently the second part of the lemma now 
follows from the first. 



6.2. Density scheme. Having established strong convergence of the velocity we 
now prove that the numerical solutions converge to a weak solution of the continuity 
equation (1.1). 

Lemma 6.4 (Convergence of the continuity approximation). The limit pair (g,u) 
constructed in (6.2) is a weak solution of the continuity equation (1.1) in the sense 
of Definition 2.4- 

Proof. The proof is essentially identical to the proof of Lemma 6.4 in [8] and is only 
included for the sake of completeness. 

Fix a test function tfi E C^°([0,T) x S7), and introduce the piecewise constant 

approximations (f> h :~ n®4>, 4>™ '■= H®(j) m , and 4> m :— ^ J tm -i 4>{t, ■) dt. 

Let us employ (/>™ as test function in the continuity scheme (3.2) and sum over 
m — 1, . . . , M. The resulting equation reads 



□ 




M 



m—1 




As in the proof of Lemma 4.6 we can rewrite this as 





(6.10) 




Lemma 4.5 tells us that 




lg h j dE (u h ■ v)-{<t> h - <t>) dS{x)dt <Chi ||V^|| 



L 2 (0,T;L 2 * (fj)) ' 



In view of Lemma 6.1, 




28 



KENNETH H. KARLSEN AND TRYGVE K. KARPER 



Summation by parts gives 

M 



•n=l Ja 

= - [ [ Q h (t-At,x)^-(tIc4>h) dxdt - [ Ql4>\dx 
J At Jn at Jn 

(j>t dxdt — / 0o0(Oj£) dx. 



h->0 



o Jn 



where (6.2), together with the strong convergence g^ go, was used to pass to 
the limit. Summarizing, letting h — ► in (6.10) delivers the desired result (2.1). □ 

6.3. Strong convergence of the density approximation. To obtain strong 
convergence of the density approximations gh, the main ingredient is a weak con- 
tinuity property of the quantity P e s(Qh, u-h). To derive this property we use (5.1) 
and a corresponding equation for the weak limit P c s{g, u). We start by deriving 
the latter. 

Let -0 € C£°(0, T) be arbitrary, set = tp(g — (g°)n) in (5.1), take the limit 
h — ► 0, and apply Lemmas 5.1 and 5.2 to obtain 

lim / / P cS (g h ,u h )ip(g - (g°)n) dxdt 

^oJo Jn t (6 n) 

lim f f (Uh) i>V x A- 1 [g - (g°) n ] dxdt. 



o Jn 



o Jn 

T 



Since the operator A 1 is self-adjoint, we can integrate by parts to obtain 

/ / d^{u h )4^ x A- 1 [q- (q°) u ] dxdt 
Jo Jn 

$ (A -1 [div x u h }) ip(g - (g°)n) dxdt 
(A-^divxH^]) ^(g-(g°}n) dxdt, 

o Jn ot 

where the last equality follows by definition of H£ (6.8). 

Next, we move ip inside the time integration and use that (g°)n is independent 
of time. This gives 

/ ^(O^A- 1 [g-(g°)n] dxdt 
o Jn 

= - f [ y+ (^A-^diVx Il c u h }) g dxdt (6.12) 

+ / / ^'(t) (A-^dW^cUh}) (g- {g°) n ) dxdt. 
Jo Jn 

At this point, we recall that (g,u) is a weak solution to the continuity equation 
(Lemma 6.4). Inserting <ft = ipA~ 1 [div x H/zUh] as test function in the weak form of 
continuity equation gives 

/ f Ju ( l 0A _1 [div a! U c u h ]) g dxdt = - [ [ ipguVj-A^ 1 [U c u h ] dxdt 
Jo Jn ut J J n 
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Setting this into (6.12) gives 

/ / ^(u/OVVxA- 1 [q- (g°)n] dxdt 

JO JQ 

= f f ^(V^A" 1 [div x U c u h ]) + V''WA- 1 [div,n C M /l ](g- (g°)n) dxdt. 
Jo Jq 

Sending h — > and applying Lemma 6.3 

lim / / d^un)^^- 1 [g- (g°) n ] dxdt 
h ^Jo Jq 

ipgu (VjA -1 [div K u}) + ip'(t) A -1 [div^ u] (g - (g°)n) dxdt. 



10 JU 

Finally, we insert this expression in (6.11) and obtain 



P e ff(<?, u)gip dxdt 
t 



o Jn 



ipgu (V^A -1 [div-ttt]) + V'(*)A _x [div x u](g - (g°)n) dxdt. 

(6.13) 



Lemma 6.5 (Effective viscous flux). Given the convergences in (6.2), 

lim / / %pP e ff(g h ,u h )g h dxdt = j j t/jP e ff(g,u)g dxdt, 
h ^°Jo Jn Jo Jn 

for allip G Cl(0,T). 

Proof. Let ip G C*(0,T) be arbitrary and set (j> = ip{gh ~ (q°)q) m (5-1) to obtain 

lim / / if/Pf (g h , u h )g h dxdt = lim / / d' t l (u h ) ^V^A" 1 [g h - (g°)n] dxdt, 
h ^°Jo Jn h ^oJo Jn 

(6.14) 

where we have also used Lemmas 5.2 and 5.1. As in (5.4) and (5.5) we can use 
summation by parts and the continuity scheme (3.2) to deduce the following equality 
for the the right-hand side: 

T r 

d h t K^V.A- 1 [g h - (q%)q] dxdt 



o Jn 

M 



= jr At f <- 1 ^(^)V,A- 1 [tf-(Q° h ) a ] +u^-V m - 1 V,A- 1 dx 
, Jn 



n 

ijjulVxA- 1 [{g h - {gl)a)] dxdt 



1 

~ At j 

M 



= J2 At i u h~ ld t v.a- 1 - (el) Q ] dx 

m=l Jn 

M . 

t At / ^'"-^r^v.A" 1 [div^™- 1 ] dx 

m=l Jn 

+ E At E / ^ lehleE W ■ - ^A-^div, w™- 1 ) dS(a 

m=l EeE h JdE 

^ J V^V, A" 1 [(^ - { e l)n)} dxdt. 
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Taking the limit h — > and applying Lemma 4.5 gives 

l-T 



lim / / (u h ) V^-V.A- 1 [g h - ( Q ° h ) Q ] dxdt 

h ^°J0 JQ 

M „ 

= lim V At / ^^uJ^VxA" 1 [div^M™- 1 ] dx 
ttV'WV.A- 1 [g-(g°)n] dxdt. 



(6.15) 



o JQ 

We will now pass to the limit in the first term on the right-hand side. 
From Lemmas 5.2, 6.1, and 6.3, we have that 

g h -» q in L°°(0, T; IS (SI)) n i 7+1 (0, T; L~< +1 (Sl)), 
u h ->u in L 2 (0,T;L 2 (fi)), (6.16) 
V X A~ X [div. T u 7l ] -> V^A" 1 [divs u\ in L 2 (0, T; L 2 {Sl)). (6.17) 

This is insufficient to pass to the limit in the desired term. However, since Uh £b 
L°°(0, T; L 2 (Sl)) n L 2 (0, T; L 2 * (0)), we can by similar arguments as in the proof 
of Lemma 6.3 deduce that V^A" 1 [div x u h ] E b L°°{0, T; L 2 {Sl)) n L 2 (0, T; i 2 * (fi)). 
Let /3 be given by 

2 1 1 
/? ~ 2 + 2*' 

Then, f3 > -j^j — e, for any e > 0. Since 2 < f3, the standard interpolation 
inequality can be applied and yields 

\f\\Un) dt ^ I* MIhvWl*-™ dt 



\Lf>(n.) Ub - / n \\J \\L 2 (n)\\J \\l 2 *(u) 

^ \\J lll ac (0,T;L 2 (O))ll/lll / 2( 0iT;L 2*( n ))- 

From this inequality, we conclude 

V.A- 1 [div, u h ] , u h e b L 4 (0, T; (6.18) 

For notational convenience, we introduce the function g^ 

g h (t, ■) = u h (t, ■) ■ V^A" 1 [div K u h (t - At, ■)] • 

Note that gh is precisely the scalar product in (6.15). From the Holder inequality 
and (6.18), we have in particular that 

g h £ b L 2 (0,T;L 2 (tt)) 

This, together with (6.16) and (6.17), tells us that 

g h -> g := itV^ A" 1 [div x u] , in L p (0, T; L p (Sl)), for any p < 2, as h -> 0. 

Hence, gngh 9Q in the sense of distributions on (0,T) x SI. This is sufficient 
to pass to the limit in the first term on the right-hand side of (6.15). By sending 
h — > in (6.15), we obtain the identity 



lim f [ d' t l (u h ) ^V.A" 1 [g h - ( e °) n ] dxdt 
h ^°Jo Jn 

= [ f tpguVxA' 1 [div x u] - u^'(t)\7 x A^ 1 [g - (g°)n] dxdt. 
Jo Jn 



(6.19) 
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Then, (6.19) in (6.14) yields 



lira / / P e s(g h ,u h )ipg h dxdt 
h ^oJo Jn 

/ / iPquVxA- 1 [div, u] - uip'itjVxA- 1 [g - (g°) n ] dxdt 
Jo Jn 

T r 

ipgu (V^A" 1 [div, u]) + ip'(t)(g - (g°)n) (A -1 div, u) dxdt 



o Jn 

T f 

/ P c s(g, u)ipg dxdt, 
la Jn 

where the last equality is (6.13). This concludes the proof. □ 

We are now in a position to prove strong convergence of the density approxima- 
tions. 

Lemma 6.6 (Strong convergence of the density). Suppose that (6.2) holds. Then, 
passing to a subsequence if necessary, 

gh — ► g a.e. in (0,T) x O. 

Proof. The proof is identical to that of Lemma 6.6 in [8] and is included for the 
sake of completeness. 

In view of Lemma 6.4, the limit (g, u) is a weak solution of the continuity equa- 
tion and hence, by Lemma 2.6, also a renormalized solution. In particular, 

(glog g) t + div, ((plog g) u) = gdiv x u in the weak sense on [0, T) x il. 

Since t h g log g is continuous with values in some Lcbcsgue space equipped 
with the weak topology, we can use this equation to obtain for any t > 

/ (glogg) (t) dx — / go log go dx = — / / g div x u dxds (6.20) 
Jn Jn Jo Jn 

Next, we specify 4>h = 1 as test function in the renormalized scheme (4.1), 
multiply by At, and sum the result over m. Making use of the convexity of zlogz, 
we infer for any m = 1, . . . , M 

/ qI log $ dx - / qI log gidx<-Y^At Qt div, < dxdt. (6.21) 
Jn Jn k=1 Jn 

In view of the convergences stated at the beginning of this section and strong 
convergence of the initial data, we can send h — > in (6.21) to obtain 

/ ( Q l°g Q ) (t) dx — / go log go dx <— / / g div, u dxds. (6.22) 
Jn ^ ' Jn Jo Jn 

Subtracting (6.20) from (6.22) gives 

/ ( £"log£? — f?logf? J (t) dx < — / / g div, u — g div, u dxds, 
Jn ^ ' Jo Jn 

for any t E (0,T). Lemma 6.5 tells us that 

rt 



, gdiv, u — gdiv x u dxds = / / g"i +1 — g^ g dxds > 0, 

o Jn M + * Jo Jn 

where the last inequality follows as in [h, 11], so the following relation holds: 



g log g — g log g a.e. in (0, T) x SI. 
Now an application of Lemma 2.1 brings the proof to an end. □ 
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6.4. Velocity scheme. 

Lemma 6.7 (Convergence of the momentum approximation). The limit triple 
(w,u,g) constructed in (6.2) is a weak solution of the velocity equation (1.2) in 
the sense of (2.3). 

Proof. Fix (v,rf) G C£°((0, T) x f2), and introduce the projections Vh = njjv, 
Vh = nf T? and v% = -fa f£_i v h dt, r% = fa J t ' m _ 1 ry^ di. 

Utilizing v™ and T7™ as test functions in the velocity scheme (3.3), multiplying 
by At, summing the result over m, and applying summation by parts, we gather 

- / / u h (t- At,x)dt (v h ) dxdt 
J At Jn 

+ \i curia, w h v h + [(/! + A) div x u h - p(gh)} div x v h dxdt = / u°vl dx, 

Jn 

/ / WhTjh - Uh curia, r\ h dxdt = 0. 
Jo Jn 

(6.23) 

In view of Lemma 2.8, Vh — ► v in L°°(0, T; W dlv,p ) for any finite p and rjh — > 
r\ in L° a (0,T;W curl ' p ). Furthermore, by Lemmas 5.2 and 6.6 p(gh) — > p(£>) in 
L a ((0, T) x 0) for any a < 7 + 1. Hence, we can send h — > in (6.23) to obtain 
that the limit constructed in (6.2) satisfies (2.2) for all test functions (v,rj) £ 
C c °°((0,T) x CI). Since C c °°((0,T) x fi) is dense in both L 2 (0, T; W curi,2 (fi)) and 
W 1,2 ^, T; L 2 {Q)) n L 2 (0, T; W div ' 2 (17)) [10] this concludes the proof. 

□ 
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